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Resumen 


López-De-La-Cruz, J., £ Francés, F. (julio-agosto, 2014). La 
variabilidad climática de baja frecuencia en la modelación no 
estacionaria de los regímenes de las crecidas en las regiones 
hidrológicas Sinaloa y Presidio San Pedro. Tecnología y 
Ciencias del Agua, 5(4), 79-101. 


El asumir estacionaridad ha sido una de las premisas 
principales en el estudio de las componentes del ciclo 
hidrológico y la piedra angular en el análisis de frecuencia 
de eventos extremos. La estacionaridad ha sido una 
hipótesis común y práctica en la planificación y gestión de 
recursos hídricos. A partir de ella se han usado métodos 
estadísticos para extraer de los datos todos los indicadores 
hidrológicos útiles para proporcionar estimaciones, donde 
estas estimaciones pueden ser refinadas año con año 
conforme los registros en las estaciones hidrométricas se 
hacen más largos. En años recientes, diversos estudios han 
demostrado que los registros hidrológicos presentan algún 
tipo de no estacionaridad, como cambios y tendencias, lo 
cual ha llevado a los expertos a un consenso, en el sentido 
de que la hipótesis de estacionaridad a nivel de cuenca se 
encuentra comprometida. Entre los principales mecanismos 
que han sido sugeridos como los inductores de cambios 
en el ciclo hidrológico de las cuencas y en la magnitud 
y frecuencia de las crecidas se encuentran los efectos de 
la intervención humana (p. ej. cambio de uso de suelo, 
embalses), el efecto de la variabilidad climática de baja 
frecuencia (p. ej. El Niño-Oscilación del Sur, Oscilación 
Decadal del Pacífico) y el cambio climático debido al 
incremento de los gases a la atmósfera. El objetivo del 
presente estudio yace en el desarrollo de un marco para 
análisis de frecuencia bajo condiciones de no estacionaridad 
por medio de los Modelos Aditivos Generalizados en 
Localización, Escala y Forma (GAMLSS por sus siglas en 
inglés). Dos diferentes aproximaciones para la modelación 
estadística no estacionaria fueron las aplicadas a los registros 
de caudales instantáneos máximos anuales en las regiones 
hidrológicas Sinaloa y Presidio San Pedro en el noroeste del 
Pacífico mexicano. Estos modelos consisten básicamente 
en el modelo con incorporación de tendencias temporales 
en los parámetros de las distribuciones paramétricas y el 
modelo con incorporación del forzamiento de la variabilidad 
climática de baja frecuencia. Los resultados en la primera 
aproximación muestran la capacidad de los modelos para 
describir la variabilidad presente en los regímenes de 
crecidas; asimismo, se observa la alta dependencia de los 
parámetros de las distribuciones paramétricas respecto 
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Stationarity has been a main premise in the study of the components 
in the hydrological cycle and the corner stone in the frequency 
analysis of extreme events. Stationarity is also a common hypothesis 
and practice in planning and managing water resources. Based on 
this assumption, statistical methods have been used to extract data 
related to all the hydrological indicators used to develop estimates. 
These estimates can be refined from year to year as the recorded 
data increase. Over recent years, a variety of studies have shown 
that hydrological records present some type of Stationarity, such 
as changes and trends, leading experts to question the stationary 
hypothesis at the basin level. The effects of human intervention (such 
as changes in land use and reservoirs), the effect of low-frequency 
climate variability (for example, El Niño-Southern Oscillation, 
Pacific Decadal Oscillation) and climate change caused by increased 
gasses in the atmosphere are among the primary mechanisms that 
have emerged as inducers of changes in the hydrological cycle of 
basins and in the magnitude and frequency of extreme floods. 
The objective of the present study is to develop a framework for 
the analysis of frequency under non-stationary conditions using 
Generalized Additive Models for Location, Scale and Shape 
(GAMLSS). Two different approaches to non-stationarity statistical 
modeling were applied to annual instantaneous peak flow records 
from the Sinaloa and Presidio San Pedro hydrological regions in the 
Pacific northwestern Mexico. These models include the model with 
temporal trends in parametric distribution parameters and the model 
with forced low-frequency climate variability. The results from the 
trend model show the ability of models to describe the variability in 
flood regimes. In addition, the parametric distribution parameters 
are observed to be highly dependent on time, which suggests a lack of 
stationarity in the flood regímes in the gauging stations studied. The 
second approach— in which the climate indices (Niño 12, Niño 3, 
Niño 3.4, SOI and PDO) that describe the behavior of low-frequency 
variability patterns are incorporated in the models as explanatory 
covariables —makes it possible to demonstrate the important role 
of the macro-scale phenomena that occur in the Pacific on the inter- 
annual variability of the flood regimes in the Pacific Mexican coast. 
In addition, a comparison of classic inference models between non- 
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del tiempo, lo cual sugiere la ausencia de estacionaridad 
en los regímenes de crecidas en las estaciones de aforo 
de estudio. En el segundo enfoque, en el cual los 
índices climáticos (Niño12, Niño3, Niño3.4, SOI y PDO) 
que describen el comportamiento de los patrones de 
variabilidad de baja frecuencia fueron incorporados como 
covariables explicativas en los modelos, permiten resaltar 
el importante papel de los fenómenos de macroescala que 
ocurren en el Pacífico, en la variabilidad interanual de los 
regímenes de las crecidas en la costa del Pacífico mexicano. 
Además, la comparación de los modelos en la inferencia 
de cuantiles entre los modelos no estacionarios respecto 
del clásico modelo estacionario muestra que las diferencias 
obtenidas asumiendo no estacionaridad y sus equivalentes 
estacionarios pueden ser importantes durante extensos 
periodos de tiempo. 


Palabras clave: análisis de frecuencia de crecidas, caudales 
máximos anuales, GAMLSS, no estacionaridad, patrones 
de variabilidad climática de baja frecuencia, Pacífico 
mexicano. 


Introducción 


La hipótesis de estacionaridad en las series 
de tiempo hidrológicas ha sido una de las 
premisas principales asumida para su estudio 
y la piedra angular en el análisis de frecuencias 
de eventos extremos (Stedinger, Vogel, éz 
Foufoula, 1993; Khaliq, Ouarda, Ondo, Gachon 
éz Bobée, 2006; Villarini et al., 2009b). De hecho, 
se puede mencionar que en la actualidad 
toda la infraestructura hidráulica con la que 
se cuenta ha sido diseñada bajo el supuesto 
de sistemas fijos; en otras palabras, que se ha 
estado trabajando en un mundo estacionario. 
La estacionaridad en una serie hidrológica se ha 
asumido en el contexto de que ésta se encuentra 
libre de tendencia o periodicidad (Salas, 1993). 
Sin embargo, diversos estudios realizados en 
las pasadas décadas han demostrado que los 
registro hidrológicos presentan algún tipo 
de no estacionaridad en forma de tendencias 
(Lettenmaier €. Wallis, 1994; Kundzewicz 
éz Robson, 2004; Stahl et al., 2010), cambios 
bruscos (Potter, 1976; Salas €: Boes, 1980) o una 
combinación de ambas (Villarini, Serinaldi, 
Smith, 8: Krajewski, 2009a). Esto ha llevado a los 
investigadores a cuestionar la hipótesis de un 
mundo estático sin influencia de forzamientos 


stationary and stationary quantiles shows that differences between 
stationarity and non-stationarity assumptions can be significant 


over long periods of time. 
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externos y posicionarse en la apremiante 
necesidad de un cambio de paradigma en este 
sentido (Khaliq et al., 2006; Milly et al., 2008; 
Sivapalan $ Samuel, 2009). Aspectos como la 
intervención humana en las cuencas, el efecto 
de la variabilidad climática de baja frecuencia, 
y el cambio climático antropogénico han sido 
sugeridos por diversos investigadores como 
los principales mecanismos de cambios en 
ciclo hidrológico y variaciones en la magnitud 
y frecuencia de las crecidas (Perreault, Haché, 
Slivitzky, € Bobée, 1999, Rasmussen, 1999; 
Zhang « Schilling, 2006; Villarini et al., 2009b). 
Dos estudios de referencia en cuanto al 
concepto de estacionaridad y no estacionaridad 
en hidrología son los presentados por Matalas 
(1997) y Koutsoyiannis (2006). 

Un trabajo reciente que colocó el tema 
de la no estacionaridad a la vanguardia de 
las ciencias hidrológicas es el publicado por 
Milly et al. (2008), en el cual enunciaron “la 
estacionaridad está muerta” y “no debe ser 
revivida”, proponiendo que la hipótesis 
de estacionaridad no debe seguir siendo 
asumida y exhortando al desarrollo de nuevos 
procedimientos estadísticos que permitan 
capturar la dinámica de evolución de las 
funciones de densidad de probabilidad sobre 
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el tiempo (Cox, Isham $ Northrop, 2002; Milly 
et al., 2008; Villarini ef al., 2009b). Ellos basan 
este urgente cambio de paradigma ante la 
necesidad de contar con la herramienta para 
adaptar el diseño y manejo de estructuras de 
control y protección a una realidad dinámica 
en las que deban adoptar nuevas estrategias y 
permitan mitigar el riesgo ante la presencia de 
eventos extremos. 

En la literatura se han propuesto diversas 
aproximaciones para la modelación de los 
regímenes de crecidas bajo condiciones de no 
estacionaridad (Olsen, Lambert, € Haimes, 
1998; Strupczewski, Singh, € Feluch et al., 
2001b; Strupczewski, Singh, $: Mitosek, 2001a; 
Cox et al., 2002; Katz, Parlange, € Naveau, 
2002; Ramesh € Davison, 2002; Cunderlink 
€ Burn, 2003; Sankarasubramanian é Lall, 
2003; Yi, Bardossy, €: Brommundt, 2006; 
Leclerc z Ouarda, 2007), donde los parámetros 
o momentos de las distribuciones de una 
distribución dada (p. ej., distribución Gumbel) 
pueden variar con el tiempo. Estos incluyen 
modelos con componentes de tendencias; 
modelos puramente estocásticos, considerando 
patrones de cambios, y modelos considerando 
covariables. Estos modelos se han planteado 
con el objetivo de no asumir las simplificaciones 
de estacionaridad en la naturaleza, permitiendo 
considerar la dinámica existente en el tiempo 
que hace de éste un sistema evidentemente más 
complejo. Unos de los esquemas abordados 
de forma reciente dentro de la modelación 
no estacionaria han sido los esfuerzos por 
aproximarse a la relación entre los extremos 
hidrometeorológicos con los patrones de 
variabilidad climática de baja frecuencia, 
lo cual se ha establecido con el objetivo de 
encontrar en los índices que describen su 
comportamiento las variables que controlen 
y mejoren los resultados en la modelación 
de eventos extremos bajo condiciones de no 
estacionaridad. Esta es una línea de estudio 
recomendada por Khalid et al. (2006), en una 
importante discusión sobre las diferentes 
aproximaciones de modelos no estacionarios. 

El objetivo principal que se plantea en el 
presente trabajo es contribuir a la modelación 


de los regímenes de crecidas bajo condiciones 
de no estacionaridad en ríos de México loca- 
lizados en el noroeste del Pacífico mexicano. 
En particular, nuestro interés radica en mo- 
delos estadísticos, incorporando tendencias 
temporales y el forzamiento de la variabilidad 
climática como covariables explicativas de los 
parámetros de las distribuciones paramétricas 
utilizadas. Los modelos no estacionarios se im- 
plementaron en el marco de los modelos GA- 
MLSS (por sus siglas en inglés), que proveen la 
flexibilidad para poder modelar los parámetros 
de las distribuciones paramétricas, incorporan- 
do covariables explicativas (Stasinopoulos éz 
Rigby, 2007). La covariables externas elegidas 
para caracterizar el efecto de la variabilidad 
climática en el estudio son índices climáticos 
que sintetizan el comportamiento de patrones 
de variabilidad climática de baja frecuencia en 
la región del Pacífico. Una vez implementadas 
las dos aproximaciones de modelos no estacio- 
narios se lleva a cabo la comparación entre el 
modelo de tendencias temporales y el modelo 
con incorporación de covariables externas res- 
pecto del clásico modelo estacionario, donde la 
comparación se basa en la estimación del cuan- 
til para una probabilidad de excedencia de 0.01 
(correspondiente al Tr de 100 años). Además, 
se lleva a cabo la predicción de crecidas para 
diferentes probabilidades de excedencia con 
base en los modelos no estacionarios ajustados, 
esto con el objetivo de observar la capacidad 
de los modelos como posibles herramientas de 
predicción. 


Efectos de la variabilidad climática 
asociada con el ENOS (El Niño- 
Oscilación del Sur) 


En las últimas décadas, un tema de amplio 
interés en la hidrología ha sido estudiar 
el impacto de la variabilidad climática en 
distintas escalas de tiempo en las componentes 
del ciclo hidrológico. El ENOS es el principal 
mecanismo modulador del clima global a 
escala de tiempo interanual. En territorio 
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mexicano, investigadores han abordado los 
efectos de la variabilidad climática sobre la 
hidroclimatología, en particular los efectos de 
las dos fases (El Niño/La Niña) del evento 
ENOS sobre los procesos hidrológicos del país. 
El Niño tiene una recurrencia promedio entre 
3 y 4 años, mientras que La Niña entre 6 y 8 
años. Estudios previos han demostrado que el 
fenómeno ENOS afecta de manera importante 
el patrón de comportamiento de las lluvias en 
México (Magaña £ Quintanar, 1997), exhibiendo 
en particular una influencia significativa en la 
evolución de las precipitaciones en el centro 
del país (Mosiño € Morales, 1988), mostrando 
no sólo repercusiones en las precipitaciones de 
verano, sino también en las invernales (Cavazos 
$ Hastenrath, 1990). En general, se puede 
mencionar que la influencia del fenómeno 
ENOS en su fase positiva (eventos El Niño) 
en la precipitación en el territorio mexicano 
es una intensificación de ésta en el noroeste 
y noreste, mientras que éstas presentan una 
disminución hacia el sur (Magaña, Pérez, 
Conde, Gat, €: Medina, 1998). Por otro lado, 
los veranos con eventos El Niño son más secos 
y cálidos que los veranos de eventos La Niña 
en casi todo el país. En cuanto a la relación de 
los fenómenos de El Niño y los regímenes de 
caudales, estudios han mostrado que durante 
los eventos El Niño se presenta un periodo 
con escurrimientos por debajo de lo normal 
en la zona del Pacífico, principalmente en el 
verano, por lo que aunque con precipitaciones 
invernales importantes, los fenómenos El Niño 
intensos han sido relacionados con sequías y 
por lo tanto con escasez de agua (Magaña et 
al., 1998). Por otro lado, en los fenómenos de 
La Niña en el Pacífico, los veranos resultan 
más húmedos, resultando esta zona una de 
las regiones donde resultan más evidentes los 
impactos de los fenómenos del ENSO. 

Ante los resultados reportados, resulta 
evidente que la función de distribución de 
probabilidad de las variables hidrológicas es 
afectada por ambas fases del ENOS, por lo que 
ante el colapso de la hipótesis de estacionaridad 
es urgente abordar la modelación estadística 


de los regímenes de crecidas desde un enfoque 
en el que se incorporen los efectos de la 
variabilidad climática (Villarini et al., 2009a). 


Caso de estudio y datos 


El presente estudio está enfocado a la 
modelación de series temporales de caudales 
instantáneos máximos anuales aforados en 
cuencas localizadas en las regiones hidrológicas 
Sinaloa y Presidio San Pedro en el noroeste 
del Pacífico mexicano (figura 1). El interés 
por abordar la modelación estadística de los 
regímenes de crecidas en sitios localizados 
en estas regiones hidrológicas yace en que en 
esta zona puede resultar más fuerte la señal 
del ENOS, esto de acuerdo con resultados 
reportados en estudios previos (Trasviña, 
Lluch, Filonov, $: Gallegos, 2004). 

Una vez revisada la base de datos en las 
dos regiones hidrológicas, se seleccionaron 
38 estaciones hidrométricas, de las cuales 25 
pertenecen a la región Sinaloa y 13 a la región 
Presidio San Pedro. Las series temporales de 
caudales máximos anuales seleccionadas para 
el estudio cuentan con al menos 30 años de 
registro (cuadro 1), las cuales fueron obtenidas 
de la base de datos de aforos BANDAS (Banco 
Nacional de Datos de Aguas Superficiales). 
Es importante mencionar que la selección de 
los eventos de caudales máximos anuales fue 
llevada a cabo considerando el año hidrológico 
del 1? de junio al 31 de mayo del año siguiente. 
En el panel superior de la figura 2 se muestra 
la distribución estacional de los eventos de 
caudales máximos anuales en las estaciones 
de aforo analizadas. En el panel superior 
izquierdo se puede observar que éstos tienden 
a concentrarse durante los meses de invierno, 
otoño y últimos meses de verano, en estaciones 
localizadas en la región Sinaloa. Por otro lado, 
un similar comportamiento en la temporalidad 
de los eventos se puede observar en las 
estaciones localizada en la región Presidio 
San Pedro (panel superior derecho), donde 
se puede distinguir una menor recurrencia 
de eventos durante los meses de invierno. En 
cuanto al régimen hidrológico en cada estación 
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Figura 1. Localización de las estaciones hidrométricas utilizadas en el estudio. 


de aforo, se evaluó con base en el coeficiente 
de caudal, en el cual se relacionó el caudal 
mensual promedio y el caudal promedio anual, 
eliminándose de esta manera el factor absoluto 
real entre los diferentes ríos. Analizando los 
gráficos del panel inferior de la figura 2, se 
puede identificar un característico régimen 
unimodal en las estaciones de las dos regiones 
con el modo en los meses de julio a octubre. 
Un régimen bimodal con un segundo modo 
atenuado en los meses de diciembre a febrero 
es el que se puede distinguir en estaciones 
localizadas en la región Sinaloa (panel 
inferior izquierdo). Estos resultados ponen de 
manifiesto el hecho de que las precipitaciones 
durante los meses de verano y otoño resultan 
ser el principal mecanismo que gobierna la 
génesis de crecidas en las estaciones de aforo 
de estudio. 

Los indicadores climáticos utilizados en 
el presente trabajo consisten en los índices 


Niñol2, Niño3, Niño3.4, Niño4, índice de 
la Oscilación del Sur (1OS) y el índice de la 
Oscilación Decadal del Pacífico (ODP). Las 
series temporales de los índices climáticos 
fueron obtenidas del Centro de Predicción 
Climática (CPC, http:/ /www.cpc.ncep.noaa. 
gov/), las cuales consisten en información a 
escala mensual. 

Estos índices proporcionan información 
acerca de la variabilidad interna de patrones 
de variabilidad climática de baja frecuencia 
en la zona del Pacífico, permitiendo así tomar 
en cuenta el impacto que puede ejercer la 
variabilidad climática de baja frecuencia en 
la magnitud y frecuencia de las crecidas. La 
selección de los índices climáticos anuales 
para su empleo como covariables explicativas 
en los modelos consistió en asumir los valores 
promedio para el periodo abril a septiembre 
previo al año hidrológico del registro de los 
eventos de crecidas. 
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Cuadro 1. Características principales de las estaciones de aforos utilizadas en el estudio. 


Caudal medio 


Núm. OR e Área de cuenca Longitud de o ua Esciciente de 
(km?) registro (años) (m/s) variación 

1 Puente Sudpacífico 10018 15 731 70 (1924-1993) 1 876.83 1.01 

2 El Bledal 10027 sy 59 (1938-1996) 283.52 0.91 

3 Canatlán 10029 2064 45 (1940-1984) 612.17 1.10 

4 Guamúchil 10031 1 645 36 (1939-1974) 593.60 1.10 

5 Palos Blancos 10033 49 (1940-1988) 1320.70 1.34 

6 Zopilote 10034 63 (1940-2002) 344.29 0.79 

7 Naranjo 10035 45 (1940-1984) 11.89 0.28 

8 Jaina 10036 61 (1942-2002) 947.12 ¡Lie 

9 Huites 10037 52 (1942-1993) 3 118.47 1.07 
10 Santa Cruz 10040 57 (1944-2000) 943.26 ¡Ly 
11 Sanalona I 10041 33 (1945-1977) 487.93 1.15 
12 Bamicori 10057 36 (1952-1987) 173.88 1.02 
13 Ixpalino 10065 48 (1953-2000) 1177.91 0.92 
14 Choix 10066 50 (1956-2005) 293.14 0.97 
15 Badiguarato 10079 41 (1960-2000) 1 003.93 1.64 
16 El Quelite 10083 33 (1961-1993) 463.26 0.94 
17 Pericos 10086 33 (1961-1993) 230.97 0.68 
18 Tamazula 10087 38 (1963-2000) 559.49 0.73 
19 Urique Il 10100 35 (1968-2002) 297.03 0.60 
20 Toahayana 10110 32 (1957-1988) 944.59 0.73 
21 Piaxtla 10111 47 (1958-2004) 1266.19 1.24 
22 Guatenipa II 10112 37 (1969-2005) 1 550.19 0.84 
23 La Huerta 10113 31 (1970-2000) 858.16 0.64 
24 Guasave 10120 31 (1975-2005) 685.10 1.03 
25 Tecusiapa 10137 47 (1958-2004) 952.70 1.08 
26 San Felipe 11008 63 (1943-2005) 126.46 LS 
27 Refugio Salcido 11010 63 (1944-2006) 122.01 1.03 
28 Peña del Águila 11011 53 (1954-2006) 91.30 1627 
29 San Pedro 11012 64 (1944-2007) 1 835.92 0.74 
30 Acaponeta 11014 62 (1946-2007) 2 087.90 1.18 
31 Baluarte II 11016 56 (1948-2003) 2 944.76 1.07 
32 Caborca 11023 56 (1950-2005) 95.30 1L2L) 
33 El Bejuco 11030 43 (1959-2001) 182.69 0.60 
34 La Ballona 11035 47 (1961-2007) 439.17 0.87 
35 El Pueblito 11036 1810 46 (1961-2006) 172.48 1.53 
36 Vicente Guerrero 11040 1 868 46 (1962-2007) 52.91 0.99 
37 Siqueiros 11058 5614 52 (1956-2007) 116.22 1.03 
38 Las Tortugas 11070 863 36 (1971-2006) 497.85 0.99 
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Figura 2. Distribución estacional de los caudales instantáneos máximos anuales (panel superior) y comportamiento 
del coeficiente de caudal mensual (panel inferior). 


Modelos aditivos generalizados en 
localización, escala y Forma (GAMLSS) 


Coincidiendo con las recientes opiniones de 
diversos autores (Milly et al., 2008; Villarini ef 
al., 2009a), al no poder garantizar la hipótesis 
de estacionaridad en el comportamiento 
de los procesos hidrológicos en las escalas 
interanuales y decadales, resulta necesario 
utilizar modelos estadísticos no estacionarios 
capaces de reproducir la variación temporal de 
los parámetros de las funciones de distribución 
probabilidad. En el presente estudio se 
utilizan los modelos GAMLSS propuestos de 
manera reciente por Rigby y Stasinopoulos 
(2005), los cuales proveen un marco flexible 
para la modelación de series temporales 
bajo condiciones de no estacionaridad. En 
los modelos GAMLSS, la respuesta de la 


variable aleatoria Y (caudales instantáneos 
máximos anuales en este trabajo) asumen una 
distribución paramétrica y sus parámetros 
pueden ser modelados como función de 
covariables seleccionadas, en este caso: tiempo 
(t) e índices climáticos (Niñol2, Niño3, 
Niño3.4, Niño4, JOS, y ODP). 

En los modelos GAMLSS, la distribución de 
la variable Y no es restringida, permitiendo el 
uso de funciones de sesgo y curtosis altos, así 
como de funciones continuas y discretas. La 
parte sistemática del modelo permite modelar 
los parámetros de localización, escala y forma 
(relacionados con la media, sesgo y curtosis) 
de la distribución de Y como lineal o no lineal, 
paramétrica o no paramétrica, a través de fun- 
ciones de suavizado (Rigby € Stasinopoulos, 
2005; Stasinopoulos éz Rigby, 2007). Para los 
modelos GAMLSS, las observaciones y, (para 
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¡ =1, 2,...n) son asumidas independientes con 
una función de distribución F,, = (y, 8,, donde 
0, = (8, e 
considerados de la distribución para localiza- 


0,,) es un vector de p parámetros 


ción, escala y forma. Usualmente p es menor o 
igual a 4, puesto que familias de hasta cuatro 
parámetros proporcionan la suficiente flexibi- 
lidad para describir los datos. Los parámetros 
de la distribución están vinculados con las 
covariables explicativas a través de g,(-) para 
k = 1,2...p, conocida como función de enlace 
monótona. En este trabajo, las funciones de 
enlace consideradas son la de identidad y la 
logarítmica. Los modelos GAMLSS envuelven 
diversos modelos, en particular nosotros 
partimos de la formulación aditiva semi- 
paramétrica, la cual está dada por: 


801) =X 181 + hy (x;) (1) 
1 


donde 6, es un vector de longitud n; 0 = (B,,,... 
P, ) es un vector de parámetros de longitud 
m; X, es una matriz de covariables de orden 
nx m; h,(x,) representa los términos aditivos 
del modelo, siendo h, una función de sua- 
vizado que evalúa el vector de covariables 
X, (Stasinopoulos € Rigby, 2007). En este 
estudio se asume que los parámetros de las 
funciones de distribución seleccionadas para el 
estudio son función del tiempo e indicadores 
climáticos. Incorporar términos de suavizado 
a los modelos (ecuación (1)), nos proporciona 
la ventaja de poder identificar no linealidades 
en la dependencia en la modelación de 
los parámetros de las distribuciones como 
función de las covariables. En este estudio, 
la dependencia suavizada está basada en las 
formulaciones matemáticas de splines cúbicos 
(Rigby € Stasinopoulos, 2005; Stasinopoulos éz 
Rigby, 2007), las cuales han mostrado buenos 
resultados en los trabajos de Villarini et al. 
(2009b) y Villarini, Smith y Napolitano (2010). 
Para el uso de las formulaciones splines cúbicos 
es necesario definir los grados de libertad, 
teniendo para el caso cuando los grados de 


libertad tienden a cero que los splines cúbicos 
tienden a una función lineal. La selección de los 
grados de libertad se lleva acabo de acuerdo con 
lo propuesto por Stasinopoulos y Rigby (2007), 
empleando los criterios de información de 
Akaike y el criterio Bayesiano de Schwarz. Por 
ello, el modelo paramétrico lineal se tiene como 
un caso particular, en el cual no se cuenta con 
términos aditivos con dependencia suavizada. 
Entonces tenemos un modelo de la forma: 


8x(01) =X Ps 0) 


donde X, B, es una combinación de estimadores 
lineales. Este modelo es incluido como el caso 
límite en la incorporación de dependencia 
suavizada. Otro caso particular es el que se 
presenta cuando todos los parámetros de las 
distribuciones resultan independientes de las 
covariables explicativas, entonces el modelo 
para los parámetros 6, se simplifica a un 
modelo estacionario g,(8,) = constante. 

En lo que respecta a la metodología para la 
selección y discriminación de las covariables 
explicativas para la descripción de los 
parámetros de las distribuciones, se siguió el 
criterio propuesto por Stasinopoulos y Rygby 
(2007). Ellos recomiendan determinar el mejor 
modelo para 6,, posteriormente para 0, hasta 
O, llevando a cabo la comparación de los 
diferentes modelos con base en los criterios 
de Akaike (CIA) y Bayesiano (CIB). Con estos 
criterios, los modelos finales proporcionan 
un balance entre exactitud y complejidad, 
tratando de evitar el sobreparametrizado de 
los modelos. Es importante mencionar que en 
ninguno de los casos, los grados de libertad 
en los splines cúbicos se tomaron mayores 
a Ln(n). Esto debido a que un aumento en la 
complejidad del modelo está vinculado con la 
extracción de información de los datos. Para 
una mayor comprensión de la teoría, ajuste y 
selección del modelo se remite a las consultas de 
Rigby y Stasinopoulos (2005), y Stasinopoulos 
y Rygby (2007). 

Una vez definida la dependencia funcional 
entre los parámetros de las distribuciones 
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y cada covariable seleccionada y los grados 
de libertad en caso de asumir dependencia 
suavizada a través de los splines cúbicos, se 
seleccionó la función de distribución F(y, 0.) 
de acuerdo con el mayor valor de la máxima 
verosimilitud. En este trabajo se seleccionaron 
cinco funciones de distribución utilizadas en 
el análisis de frecuencias de crecidas (cuadro 
2): Gumbel (GU), Lognormal (LNO), Weibull 
(WED), Gamma (GA) y Gamma Generalizada 
(GG), de las cuales las cuatro primeras son de 
dos parámetros y la última de tres (Stedinger et 
al., 1993; Escalante $ Reyes, 2002). 

En la ausencia de estadísticos para evaluar la 
calidad del ajuste de los modelos seleccionados 
en su conjunto, la verificación se realizó de 
conformidad con las recomendaciones de Rigby 


y Stasinopoulos (2005) mediante el análisis de la 
normalidad e independencia de los residuales 
de cada modelo. Para ello, se examinaron 
los cuatro primeros momentos estadísticos 
de los residuales, el autocorrelograma y los 
coeficientes de correlación de Filliben, así como 
la inspección visual de los gráficos residuales ús. 
respuesta, q-q plot y q-q plot sin tendencia. Esta 
acción garantiza que los modelos seleccionados 
pueden explorar adecuadamente la parte 
sistemática, siendo la información restante 
(residual) ruido blanco (señal aleatoria). Todo 
los cálculos se llevan a cabo en el paquete 
GAMLSS (Rigby é Stasinopoulos, 2005; 
Stasinopoulos € Rigby, 2007), implementados 
en la plataforma R (R Development Core Team, 
2008). 


Cuadro 2. Resumen de las funciones de densidad de probabilidad consideradas en la modelación de las series 


de crecidas máximas anuales. 


Función de densidad de probabilidad 
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Resultados 
Implementación de los modelos en GAMLSS 


En esta sección se presentan y analizan los 
resultados obtenidos en la modelación de los 
regímenes de crecidas en las estaciones de 
estudio bajo condiciones de no estacionaridad. 
Se implementaron tres diferentes tipos de 
modelos: a) Modelo de tendencia, b) Modelo 
de covariables y c) Modelo estacionario. 


Modelación de los regímenes de crecidas con 
incorporación de tendencias temporales: 
modelo de tendencia 


La primera aproximación en la modelación 
de los regímenes de crecidas bajo condiciones 
de no  estacionaridad consistió en la 
implementación de modelos en los cuales se 
incorporó el tiempo como única covariable 
explicativa. En el cuadro 3 se resumen los 
modelos de tendencias temporales ajustados 
a las 38 series temporales de crecidas de los 
sitios de estudio, por medio de los cuales se 
intentó describir la variabilidad interanual 
presente en los regímenes de crecidas en el 
tiempo. El ajuste de los modelos no presenta 
una distribución como la descriptora del total 
de las series temporales, sólo es de destacar 
que las distribuciones Gamma y Lognormal 
son las que mejor describen el comportamiento 
de los regímenes de crecidas de acuerdo con 
las distribuciones utilizadas en el estudio. 
En cuanto a la dependencia, así como al tipo 
de dependencia de los parámetros de las 
distribuciones en el tiempo, se observa que el 
parámetro 6, presenta una fuerte dependencia 
respecto del tiempo, encontrando dicha de- 
pendencia en la mayoría de los modelos de 
tipo suavizada. Una fuerte dependencia es 
la observada también con el parámetro 0»; 
no obstante, resulta evidente que ésta resulta 
más atenuada. Es importante destacar que 
se pueden distinguir tanto dependencia de 
tipo suavizada como dependencia lineal. En 
general, es posible mencionar que los sitios 


de estudio presentan una fuerte dependencia 
del tiempo (34 estaciones), donde 15 sitios 
exhibieron dependencia en los dos parámetros, 
mientras que sólo en cuatro sitios se detectó 
una independencia total de los parámetros. 

La modelación de las series temporales 
de caudales máximos anuales empleando el 
tiempo como covariable explicativa de los pa- 
rámetros de las distribuciones paramétricas 
muestra la habilidad de los modelos ajustados 
de capturar la variabilidad exhibida en los re- 
gímenes de crecidas. En la figura 3 se muestran 
los resultados de la modelación de los cauda- 
les máximos anuales con la aproximación del 
modelo no estacionario de tendencia temporal 
en seis sitios representativos de los 38 analiza- 
dos. Analizando los resultados expuestos en la 
figura 3, se puede observar una cierta similitud 
en la variación en el tiempo de los cuantiles es- 
timados para diferentes probabilidades de no 
excedencia en las estaciones que se encuentran 
más cercanas a la costa, 10027, 10034, 11012 y 
11058; en ellas se aprecia un patrón creciente en 
la magnitud de las crecidas durante el periodo 
1940-1970, después se identifica una tendencia 
decreciente; no obstante, no es un patrón ge- 
neralizado, ya que se puede observar que en 
algunas estaciones comienza en la década de 
los ochenta (10034 y 11012), mientras que en 
otros casos se prolonga hasta la década de los 
setenta (10027 y 10011058). En cuanto al com- 
portamiento en estaciones localizadas en la 
parte más interior de las regiones de estudio, 
se observan resultados un tanto diferentes a 
los observados en las estaciones anteriores. 
Analizando las estaciones presentadas en la fi- 
gura 3 (10066 y 11008), se puede ver en estas 
estaciones la presencia de una tendencia cre- 
ciente en la magnitud de las crecidas posterior 
a la década de los cincuenta, mientras que la 
estación 10066 exhibe un decremento en el pe- 
riodo comprendido entre los últimos años de 
los ochenta hasta finales de los años noventa. 
Un aspecto a destacar en los resultados obte- 
nidos en la implementación de los modelos de 
tendencia yace en que la variabilidad temporal 
en los cuantiles resulta más importante en los 
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Cuadro 3. Resumen de los modelos de tendencia ajustados y el tipo de dependencia entre el tiempo y los parámetros de las 
distribuciones, así como de lo estadístico de los residuales de los modelos y el coeficiente de correlación de Filliben: 
cs(-) indica que la dependencia es modelada por medio de los splines cúbicos, “t” indica dependencia lineal y “--““se refiere 
a que el parámetro es independiente del tiempo (estacionario). 


Estación Distribución 8, 8, Media Varianza peris AIC SBC 

10018 LNO cs(t) cs(t) 0.0163 1.0142 0.9914 1113.74 1 136.22 

10027 LNO cs(t) cs(t) 0.0009 1.0172 019955 782.77 803.55 

10029 LNO - cs(t) -0.0551 1.0201 0.9948 655.49 666.21 

10031 GA - cs(t) -0.0146 1.0053 0.9807 518.49 527.82 

10033 LNO cs(t) cs(t) -0.0036 1.0203 0.9946 771.51 790.63 

10034 WEI cs(t) - -0.0014 1.0187 0.9969 862.75 875.62 

10035 LNO cs(t) cs(t) 0.0137 1.0225 0.9855 188.83 206.89 

10036 LNO P - 5.518e-17 1.0166 0.9794 O2D:D2 931.86 

10037 LNO cs(t) cs(t) 0.0002 1.0196 0.9493 933.04 952.56 

10040 LNO cs(t) cs(t) 0.0197 1.0174 0.9892 838.33 858.76 

10041 LNO cs(t) - -1.660e-14 1.0312 0.9943 463.31 472.29 

10057 GA cs(t) - -0.0116 1.0198 0.9732 433.64 443.14 

10065 GA cs(t) - 0.0005 1.0220 0.9870 766.50 777.73 

10066 LNO cs(t) cs(t) -0.0009 1.0204 0.9812 642.90 662.03 

10079 LNO t - 2.397e-16 1.0250 0.9717 635.09 640.24 

10083 GA cs(t) cs(t) 0.0076 1.0362 0.9899 483.25 498.22 

10086 GA t - 0.0062 1.0347 0.9865 416.93 421.42 

10087 WEI cs(t) cs(t) 0.0403 1.0897 0.9772 542.74 pO9ila 

10100 WEI - t -0.0058 1.0855 0.9828 461.64 466.31 

10110 WEI - - -0.0003 1.0343 0.9833 501.42 504.35 

10111 LNO - - 1.078e-16 1.0217 0.9879 755.92 759.62 

10112 LNO p - -4.548e-17 1.0277 0.9691 614.64 619.47 

10113 GA cs(t) - -0.0031 1.0318 0.9892 468.29 476.89 

10120 GA cs(t) cs(t) -0.0083 1.0116 0.9887 461.37 475.72 

10137 GA cs(t) cs(t) -0.0031 1.0174 0.9889 741.31 759.81 

11008 GA - l -0.0092 1.0058 0.9819 720.97 727.45 

11010 GA cs(t) cs(t) -0.0151 0.9930 0.9952 730.97 752.40 

11011 WEI - - -0.0067 1.0055 0.9526 566.32 570.26 

11012 LNO cs(t) - 1.760e-14 1.0158 0.9836 1 059.72 1 072.67 = 

11014 LNO cs(t) cs(t) 0.0132 1.0162 0.9893 1 045.79 1 067.06 5 

11016 GA - cs(t) -0.0376 0.9684 0.9467 1 006.16 1 018.31 2 

11023 WEI - t -0.0124 1.0100 0.9864 615.72 621.79 E 

11030 WEI pa t 0.0079 0.9806 0.9816 512.17 519.22 3 

11035 GA cs(t) cs(t) -0.0041 1.0230 0.9822 631.01 649.08 y 

11036 LNO - - 1.461e-17 1.0222 0.9496 538.57 542.23 E 

11040 GA - cs(t) -0.0393 0.9728 0.9805 456.79 467.77 7 

11058 LNO cs(t) - -1.563e-14 1.0196 0.9914 830.82 842.53 S 

11070 GA cs(t) - 0.0211 1.0439 0.9519 497.99 507.49 a 
3 
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Figura 3. Resumen de los resultados obtenidos con el modelo de tendencia en seis estaciones representativas. Los resultados 


muestran los cuantiles estimados para diferentes probabilidades de no excedencia (0.025, 0.25, 0.50, 0.75, 0.95). 


cuantiles más altos (0.75 y 0.99), mientras que 
ésta resulta más atenuada en la mediana (cuan- 
til de 0.5). 

Como forma de evaluar la eficacia en el 
ajuste de los modelos implementados en esta 
aproximación, se llevó a cabo la verificación 
de la normalidad de los residuales para cada 
modelo. En el cuadro 3 se presenta un resumen 
de los momentos estadístico de los residuales, 
así como la estimación del coeficiente de 
correlación de Filliben entre los cuantiles, los 
cuales no muestran en general desviaciones 
importantes de la normalidad. Además de 
la revisión de los momentos estadísticos 
de los residuales, se verificaron los gráficos 
de cuantiles q-q plot sin tendencia, en los 
cuales se muestran las desviaciones de las 
observaciones de la muestra respecto a la recta 
normal de forma que si la muestra procede 
de una población normal, los puntos deben 
fluctuar alrededor de 0 y sin seguir un patrón 
determinado, quedando dentro de los límites 
de confianza (95% de nivel de confianza); en 
caso contrario, se presentaría la ausencia de 
normalidad de los residuales. 

Dentro de los puntos a destacar en la im- 
plementación de este primer modelo no es- 


tacionario se tiene la presencia de una fuerte 
dependencia de los parámetros de las distribu- 
ciones seleccionadas en relación con el tiempo. 
Aunado a lo anterior, es de destacar también la 
presencia de no linealidad en la relación. Esto 
se establece ante la necesidad de incorporar las 
formulaciones matemáticas de suavizado para 
modelar la dependencia de los parámetros en 
función del tiempo. Es importante mencionar 
que en el ajuste de los modelos, asumiendo 
dependencia suavizada en los parámetros, se 
determinaron los grados de libertad, respetan- 
do el principio de parsimonia; para ello se re- 
dujeron los mismos, en caso de no afectar en 
forma significativa la eficacia en el ajuste de los 
modelos. 

La modelación de los caudales máximos 
anuales bajo condiciones de no estacionaridad 
en este primer escenario planteado nos ha 
permitido la reconstrucción de los cuantiles 
para diferentes probabilidades de no exced- 
encia, modelando los parámetros de las 
distribuciones paramétricas, empleando el 
tiempo como covariable explicativa. El análisis 
de los gráficos permite apreciar que en general 
los modelos ajustados describen de modo 
adecuado la variabilidad que se observa en las 
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series de caudales máximos anuales. Por otro 
lado, con estos modelos es patente que un alto 
porcentaje de las series temporales exhiben 
desviación del supuesto de estacionaridad, 
siendo evidente que tanto la media como la 
varianza de las distribuciones se ven afectadas 
por los cambios de los parámetros en el tiempo. 


Modelación de los regímenes de crecidas, 
incorporando los efectos de la variabilidad 
climática: modelo de covariables 


La implementación de los modelos de tendencia 
temporal en el apartado anterior exhibe de 
manera evidente las desviaciones del supuesto 
de estacionaridad en el comportamiento de los 
regímenes de crecidas en las estaciones de aforo 
de estudio. En este apartado se implementan 
modelos no estacionarios con la incorporación 
de los índices climáticos como covariables 
explicativas de los parámetros. Estos modelos 
se implementan con el objetivo de explorar la 
factibilidad de describir de forma adecuada 
las no estacionaridades observadas en los 
modelos de tendencia temporal a través de los 
indicadores climáticos seleccionados en este 
estudio. 

En los últimos años, uno de los temas en 
que se han enfocado diversos investigadores 
ha sido estudiar la influencia que ejerce el 
fenómeno ENOS (El Niño-Oscilación del Sur) 
en los procesos hidrológicos, el cual representa 
el fenómeno océano-atmósfera dominante en 
el Pacífico tropical (Cane, 1992). Los resultados 
expuestos en estos estudios (Amarasekera, 
Lee, Williams é€z Eltahir, 1997; Jain é€z Lall, 
2001; Poveda et al., 2002; Zhang, Xu, Jiang, 
£ Wu, 2007) han puesto en evidencia que 
las fases extremas de este fenómeno pueden 
estar fuertemente vinculadas con los mayores 
episodios de crecidas y sequías en diversas 
partes del mundo. 

El cuadro 4 presenta un resumen de los 
modelos obtenidos, incorporando los efectos 
de la variabilidad climática de baja frecuencia; 
en él se muestran las funciones de distribución 
de mejor ajuste, así como las covariables 
significativas para cada parámetro. En general, 


los resultados en la selección de los modelos 
exhiben la significancia estadística de losíndices 
climáticos como covariables explicativas en la 
evolución de los regímenes de crecidas en los 
sitios de estudio. En particular, es de destacar 
la mayor significancia de los índices Niño3.4, 
Niño4 y ODP como covariables entre los seis 
índices utilizados. Por otro lado, en cuanto a la 
dependencia que se presenta en estos modelos, 
se puede observar en el cuadro 4 que el 
parámetro 0, presenta una fuerte dependencia 
respecto de los índices climáticos, mientras 
que una más débil dependencia es la que se 
puede observar con el parámetro 6,, donde 
18 sitios presentaron independencia de los 
índices climáticos. En cuanto a los resultados 
observados con el parámetro 0, en aquellos 
modelos donde la distribución de mejor ajuste 
resultó la distribución GG, éste presentó una 
total independencia, lo cual podría mostrar el 
hecho de una nula influencia de la variabilidad 
climática en los momentos de mayor orden. 
En cuanto al tipo de dependencia, se puede 
observar en el cuadro 4 que no se distingue un 
patrón generalizado, siendo sólo de destacar 
la necesidad de incorporar dependencia a 
través de las formulaciones matemáticas de 
suavizado. Este es un punto importante del 
análisis, ya que resalta la necesidad de asumir 
no linealidad en la modelación de los regímenes 
de crecidas, incorporando los índices climáticos 
como covariables. 

La figura 4 resume los resultados en la im- 
plementación de los modelos no estacionarios 
con la incorporación de covariables externas 
(adicionales al tiempo) en seis estaciones 
representativas. Se puede observar que los 
modelos tienden a reproducir de modo 
adecuado la variabilidad interanual de los 
regímenes de crecidas; de hecho, se puede 
destacar la mejor descripción respecto de los 
resultados observados con los modelos de 
tendencia temporal. En cuanto a la revisión 
de la normalidad de los residuales para el 
adecuado ajuste de los modelos, se puede 
mencionar que se ajustan razonablemente bien 
de acuerdo con lo observado en los momentos 
estadísticos de los residuales, así como en los 
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Cuadro 4. Resumen de los modelos con covariables externas y el tipo de dependencia entre los índices climáticos y los 


parámetros de las distribuciones, así como del coeficiente de correlación de Filliben para los residuales de los modelos: cs() 


indica que la dependencia es modelada por medio de los splines cúbicos, otro indica dependencia lineal y 


parámetro es independiente de los índices climáticos (estacionario). 


“24m 


se refiere a que el 


Estación | Distribución 8, 8, 8, ri 
10018 LNO Niño3 + cs(Niño3.4) + Niño4 - 0.9876 
10027 NO cs(Niño3) + Niño3.4 IOS 0.9950 
10029 LNO cs(Niño 12) + Niño 3 + Niño4 - 0.9859 
10031 GA cs(ODP) == 0.9949 
10033 LNO - ODP 0.9788 
10034 GG Niño3.4 + cs(Niño4) Niño4 - 0.9915 
10035 LNO Niño3.4 + Niño4 + IOS + ODP - 0.9747 
10036 LNO - ODP 0.9896 
10037 LNO cs(Niño3.4) + Niño4 - 0.9647 
10040 LNO Niño3 + cs(ODP) - 0.9774 
10041 LNO cs(Niño3.4) + cs(ODP) Niño3 0.9891 
10057 GA cs(Niño3.4) ODP 0.9921 
10065 GG - Niño4 + IOS - 0.9907 
10066 LNO cs(Niño3.4) + cs(Niño4) -- 0.9818 
10079 LNO - - 0.9768 
10083 GA - -- 0.9829 
10086 GA cs(Niño4) - 0.9833 
10087 WEI Niño4 + IOS + ODP cs(Niño3.4) + ODP 0.9873 
10100 GG cs(Niño4) + IOS + ODP cs(Niño3) + Niño3.4 -- 0.9793 
10110 WEI - - 0.9833 
1011 LNO cs(Niño3.4) cs(SOD) 0.9913 
10112 LNO - - 0.9832 
10113 GA - - 0.9856 
10120 GG cs(Niño3.4) + Niño4 ODP = 0.9963 
10137 GA - - 0.9727 
11008 GA Niñol12 + cs(Niño4) es(ODP) 0.9953 
11010 GG cs(Niño3.4) + IOS - - 0.9943 
11011 WEI cs(Niño12) == 0.9677 
11012 LNO cs(Niño3.4) cs(Niño12) + cs(IOS) 0.9894 
11014 LNO Niño3 + cs(Niño3.4) cs(Niño3.4) 0.9848 
11016 GA cs(Niño3.4) Niño4 + IOS 0.9808 
11023 GA Niño3.4 + cs(1OS) = 0.9875 
11030 WEI ODP Niño12 + cs(Niño3) + Niño4 0.9839 
11035 GA Niño12 + Niño3 + ODP cs(Niño3) 0.9834 
11036 LNO Niño12 + cs(ODP) - 0.9839 
11040 GG Niño3 ODP == 0.9957 
11058 LNO cs(Niño3.4) + Niño4 Niño3.4 0.9887 
11070 GA Niño3.4 + cs(Niño4) Niño4 + cs(ODP) 0.9786 
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Figura 4. Resumen de los resultados obtenidos con el modelo de covariables externas en seis estaciones representativas. 


Los resultados muestran los cuantiles estimados para diferentes probabilidades de no excedencia (0.025, 0.25, 0.50, 0.75, 0.95). 


gráficos de cuantiles. En la figura 5 se presentan 
los gráficos de cuantiles q-q plot sin tendencia 
obtenidos en la revisión de la normalidad de 
los residuales. 

Retomando la inspección de los resultados 
en la figura 4, es posible observar que la 
estimación de la mediana (cuantil para una 
probabilidad de excedencia 0.50) tiende a 
presentar un menor efecto de la variabilidad 
climática; sin embargo, los efectos son 
importantes en la estimación para los cuantiles 
más altos. Otro punto a destacar es que a 
pesar del adecuado ajuste de los modelos, se 
puede observar que existe un cierto ruido en 
los resultados del modelo, lo cual refleja la no 
linealidad en la respuesta de los regímenes 
de crecidas al forzamiento de la variabilidad 


climática natural. 


Comparación entre los modelos no 
estacionarios y el modelo estacionario 


El estudio de las crecidas en la hidrología 
operacional tiene como objetivo principal 
estimar eventos de crecida para una probabi- 
lidad de excedencia definida a priori para 
poder obtener mapas de inundación, diseñar 


medidas de protección o establecer planes de 
gestión del riesgo de inundación. De hecho, 
en México y en diversas partes del mundo, las 
legislaciones sobre el riesgo de inundación se 
basan en el análisis de frecuencia de crecidas 
para estimar crecidas de diseño asociadas con 
diversos periodos de retorno (p. ej., 20, 50 y 100 
años), donde dichos periodos de retorno están 
vinculados con la necesidad de seguridad en la 
estructura. 

La figura 6 muestra los resultados del aná- 
lisis de frecuencia de crecidas en condiciones 
estacionarias y bajo condiciones de no estacio- 
naridad para una probabilidad de excedencia 
de 0.01 (esto es, un periodo de retorno de 100 
años). La inspección de los gráficos resalta el 
problema de asumir estacionaridad en la es- 
timación de los eventos de crecida. Se puede 
observar que los modelos no estacionarios in- 
dican la existencia de importantes periodos en 
los cuales la crecida estimada bajo condiciones 
de no estacionaridad está por encima de la ob- 
tenida, asumiendo el modelo estacionario. En- 
focándonos en el análisis de los resultados para 
la estación 10027, se puede ver que asumien- 
do el modelo de tendencia, la crecida para una 
probabilidad anual de excedencia de 0.01 du- 
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Figura 5. Gráficos de cuantiles q-q plot sin tendencia de los residuales de los modelos con covariables externas 


en las estaciones representativas analizada. Las dos líneas punteadas corresponden a los límites de confianza para el 95%. 
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Figura 6. Cuantiles estimados de las crecidas máximas anuales con 0.01 probabilidad de excedencia para las seis estaciones 


representativas, basados en los tres diferentes modelos implementados. 


rante los 59 años de registro ha oscilado desde tras que asumiendo el modelo de covariables, 
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el mínimo valor es de 295.51 m*/s en 1988 y el 
máximo de 3 508.08 m*/s en 1969. 


un valor mínimo de 640.40 m*/s en 19% a un 
valor máximo de 1 730.52 m*/s en 1966; mien- 
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Analizando los resultados en la modelación 
de los eventos de crecidas en el escenario no 
estacionario con el modelo de covariables 
y el escenario estacionario para la estación 
10027 que se muestran en la figura 7, es 
posible observar que la mediana tiende a ser 
subestimada en el modelo estacionario para 
valores altos negativos del índice Niño3.4, 
mientras que el modelo sobrestima la mediana 
para valores altos positivos. Estos valores son 
más significativos para el percentil del 95%, 
como se puede ver en la figura 7. Es importante 
mencionar que se aprecia una compleja relación 
entre el índice climático y los valores de caudales 
máximos anuales, lo cual puede ser reflejo de 
que solamente una parte de la variabilidad de 
los regímenes de crecidas es explicada por los 
patrones de baja frecuencia. Esto, sin duda, 
lleva a la necesidad de considerar los efectos 
que tienen factores de origen atropogénico en 
los regímenes de crecidas, ya sea a través de 
cambios en las cuencas (p. ej. cambios de uso 
de suelo) o directamente en el ciclo hidrológico 
por medio de la construcción de embalses. 

En general, los resultados presentados en 
este apartado demuestran que las inferencias 
de los eventos de crecidas bajo condiciones de 


no estacionaridad pueden experimentar cam- 
bios importantes. Estos resultados refuerzan 
el cuestionamiento reciente de la hipótesis de 
estacionaridad, y nos llevan a plantearnos el 
posible colapso de la hipótesis de estacionari- 
dad en el estudio de las crecidas. 


Los conceptos de periodo de retorno y riesgo 
en el contexto no estacionario 


En general, los resultados obtenidos en la 
modelación de los regímenes de crecidas 
bajo condiciones no estacionarias nos llevan 
a plantear la necesidad de utilizar modelos 
alternativos para asumir la dinámica de 
la naturaleza en vez del clásico análisis de 
frecuencias de crecidas. Sin embargo, en un 
mundo no estacionario resulta necesario 
redefinir los conceptos de periodo de retorno 
y riesgo, lo anterior como consecuencia de 
los cambios que experimentará la función 
de densidad de probabilidad en el tiempo 
(Sivapalan €: Samuel, 2009; Salas éz Obeysekera, 
2013). 

La figura 8 muestra la estimación del 
periodo de retorno bajo condiciones de no 
estacionaridad con base en los modelos 
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Figura 7. Estimación de la mediana y los percentiles para el 5 y 95% para el modelo estacionario (línea negra) 


y el modelo de covariables (línea gris) para la estación 10027, donde se presentan graficados frente al índice Niño 3.4. 
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estadísticos no estacionarios ajustados para 
la serie temporal de crecidas en la estación 
10027. Se pude observar que la probabilidad de 
excedencia cambia año con año, ante lo cual es 
necesario articular un marco donde se puedan 
extender los conceptos de periodo de retorno 
y riesgo, los cuales, en su definición actual, 
pierden sentido. 


Evaluación del modelo no estacionario de 
covariables en predicción 


En este apartado se evalúa la faceta de los 
modelos de covariables como herramientas 
predictivas. Sólo considerar el modelo no 
estacionario con covariables yace en que los 
modelos de tendencia no pueden describir de 
forma adecuada el comportamiento a futuro 
de los eventos de crecidas, lo anterior debido a 
que las tendencias temporales pueden cambiar 
en el corto o largo plazos. La predicción se 
realizó a partir de la fecha última de registro 
para cada una de las estaciones con base en los 
modelos ajustados. 

La figura 9 muestra los resultados obteni- 
dos en el empleo de los modelos con covaria- 
bles como herramientas de predicción de las 
crecidas para diferentes probabilidades de ex- 
cedencia. Como se puede observar, los índices 
climáticos que caracterizan el comportamiento 


de los patrones de baja frecuencia pueden ser 
de gran ayuda para capturar los cambios en la 
frecuencia de crecidas en las estaciones de es- 
tudio. Sin duda, los modelos con incorporación 
de covariables adicionales al tiempo resultan 
una opción más atractiva en relación con los 
modelos con la incorporación de tendencias 
temporales. Esto como consecuencia de que los 
modelos de tendencia pueden ignorar cambios 
en la frecuencia de las crecidas en el periodo de 
predicción. 

Es importante tomar en cuenta que los 
índices utilizados para caracterizar el impacto 
de la variabilidad climática tienen limitaciones 
en la representación de procesos físicos 
adicionales, que pueden afectar la frecuencia 
de crecidas, por lo que a pesar de su potencial 
utilidad, los valores obtenidos en predicción 
deben tomarse con mucha precaución. Por otro 
lado, se debe considerar el hecho de que los 
índices que describen el comportamiento de 
los patrones de variabilidad de baja frecuencia 
son una simplificación de la circulación general 
de la atmósfera. 


Discusión 
En este estudio se ha examinado la no esta- 


cionaridad en la distribución de los caudales 
máximos anuales registrados en 38 estaciones 
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Figura 8. Estimación del periodo de retorno bajo condiciones de no estacionaridad con el modelo de tendencias 


(panel izquierdo) y el modelo de covariables (panel derecho) para la estación 10027. 
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Figura 9. Resultados de la modelación de las crecidas máximas anuales en las estaciones 10027, 10066, 1008 y 11058 con el 


modelo de covariables. Los resultados en el periodo de ajuste de los modelos se muestran en líneas continuas, mientras 


en el periodo de predicción con líneas continuas. 


de aforo localizadas en las regiones hidroló- 
gicas Sinaloa y Presidio San Pedro (Pacífico 
mexicano). Se implementaron dos aproxima- 
ciones de modelos bajo condiciones de no es- 
tacionaridad: el modelo de tendencia temporal 
y el modelo de covariables. Los modelos no es- 
tacionarios se desarrollaron en el marco de los 
modelos GAMLSS, los cuales mostraron la flexi- 
bilidad para la modelación probabilística de las 
series temporales de caudales máximos anua- 
les en un contexto no estacionario, así como la 
capacidad para modelar la dependencia de los 
parámetros de las distribuciones respecto de 
covariables explicativas. La implementación 
del modelo de tendencia temporal en el cual el 
tiempo fungió como única covariable permitió 
identificar que los regímenes de crecidas en un 
importante número de las estaciones de estu- 
dio experimentan importantes desviaciones 
del supuesto de estacionaridad. Analizando el 


tipo de dependencia observado en el modela- 
do de los parámetros como función del tiempo 
se puede mencionar que se presentaron tanto 
casos de dependencia lineal como dependen- 
cia no lineal a través de las formulaciones de 
suavizado. 

La aplicación de los modelos no estacio- 
narios, incorporando los efectos de la variabi- 
lidad climática, exhibió la significativa influen- 
cia que puede ejercer el fenómeno ENOS en 
la variabilidad interanual de los regímenes de 
crecidas. En particular, es de destacar que los 
altos valores de crecidas en las estaciones de 
estudio están vinculados con la fase negativa 
del ENOS (eventos La Niña), mientras que 
los valores bajos de crecidas lo tienden a estar 
con la fase positiva (eventos El Niño). Los 
resultados mostraron que incorporando de 
manera simple los efectos de la variabilidad 
climática a través de los índices climáticos, 
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éstos pueden describir de manera razonable 
las no estacionaridades observadas en el 
comportamiento de los regímenes de crecidas. 

Las aproximaciones de modelos no estacio- 
narios mostraron que las tendencias y los efectos 
de la variabilidad climática afectan tanto la 
media como la varianza de las distribuciones. 
Un punto importante a destacar es que se 
observó que aquellos modelos en los cuales 
se involucran las formulaciones de suavizado 
(splines cúbicos) tienden a reproducir mejor 
la dispersión de las crecidas, al presentar 
una mayor flexibilidad. Sin embargo, estos 
tipos de modelos que proporcionan un buen 
ajuste y flexibilidad son altamente sensibles 
a los cambios en la evolución de las variables 
predictivas, por lo que deben ser empleados 
con precaución. Otro punto a destacar en el 
uso de las formulaciones de suavizado es que 
el criterio de Akaike tiende a proporcionar 
modelos muy complejos, mientras que el 
criterio Bayesiano, modelos más simples. 
Por ello, es importante en la optimización de 
los grados de libertad en los modelos, seguir 
un procedimiento en el cual se cumpla con el 
principio de parsimonia, ya que el aumento en 
la complejidad del modelo está vinculado con 
la extracción de información de los datos. 

Los resultados del análisis de frecuencia 
de crecidas con los modelos no estacionarios 
mostraron que para un evento de crecida con 
una probabilidad anual de excedencia de 0.01 
(correspondiente a un periodo de 100 años), 
las variaciones obtenidas respecto del clásico 
modelo estacionario son muy importantes, con 
extensos periodos en los cuales el valor estima- 
do bajo condiciones no estacionarias es mayor 
que el estimado asumiendo estacionaridad. 
Estos resultados tienen un importante efecto en 
la práctica hidrológica y son una evidencia de 
que la simplificación de estacionaridad puede 
llevar a asumir un riesgo mayor al planificado 
en el diseño de estructuras hidráulicas. 

La implementación de los modelos no 
estacionarios como herramientas de predicción 
muestra que los modelos con covariables 


se presentan como una mejor alternativa 
respecto de los modelos de tendencias, ya que 
las tendencias a corto y largo plazos pueden 
cambiar por el efecto de la variabilidad 
climática y actividades humanas. Sin embargo, 
es importante considerar que en los modelos 
de covariables existen factores en el proceso de 
producción de escorrentía que no son tomados 
en cuenta. Por ello, la consolidación de los 
modelos como herramienta de predicción 
requiere de trabajos más exhaustivos. 

Se hace necesario llevar a cabo un análisis 
de varianza entre los índices climáticos, esto 
ante la alta multicolinealidad que puede existir 
entre los indicadores climáticos empleados, 
dado que describen un mismo fenómeno, 
lo cual permitirá obtener modelos que se 
ajusten al principio de parsimonia. Por otro 
lado, es importante mencionar que no se 
consideran en este estudio factores de origen 
antropogénico (p. ej., cambios de uso de suelo, 
construcción de embalses) que pueden afectar 
los regímenes de crecidas. Deben abordarse 
estudios más exhaustivos en los que se 
incorporen dichos efectos. Finalmente, aunque 
de las distribuciones paramétricas utilizadas, 
la distribución Gamma y Lognormal parecen 
proporcionar resultados aceptables, resulta 
necesario incorporar en estudios futuros 
distribuciones adicionales, como la General 
de Valores Extremos o distribuciones de dos 
poblaciones como la TCVE (Two Components 
Extreme Value). 


Conclusiones finales 


La hipótesis de estacionaridad en el estudio de 
las crecidas se encuentra comprometida y es 
necesario implementar métodos estadísticos 
alternativos que permitan asumir que las 
funciones de distribución de probabilidad se 
ven afectadas por el impacto de factores de 
origen antropogénico y climático. EL ENOS 
no es causado por el cambio climático; sin 
embargo, una pregunta de gran calado es si los 
eventos del ENOS en un escenario de cambio 
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climático serán más intensos y frecuentes. A la 
luz de los resultados, una tarea fundamental 
es incorporar los efectos de factores de 
origen antropogénico ante las evidencias 
que muestran que los impactos de estos 
factores se han intensificado en las últimas 
décadas. Aunque resultó evidente el potencial 
de los índices climáticos como covariables 
explicativas de la variabilidad en los regímenes 
de crecidas en las estaciones de estudio, es 
necesario que estas investigaciones vayan de 
la mano con las mejoras que están llevando 
a cabo los meteorólogos en el pronóstico de 
los índices, a fin de consolidarlos como un 
instrumento útil en la predicción a mediano 
y largo plazos. Además de lo anterior, se debe 
considerar el hecho de que los índices que 
describen el comportamiento de los patrones 
de variabilidad de baja frecuencia son una 
simplificación de la circulación general de la 
atmósfera, por lo que los valores en predicción 
deben tomarse con precaución. Por último, 
aunque los resultados demuestran la influencia 
del ENOS en la frecuencia e intensidad de los 
regímenes de crecidas en el noroeste del Pacífico 
mexicano, estos resultados no son extrapolables 
a otros sitios del territorio mexicano, ya que 
de acuerdo con los resultados reportados en 
estudios previos sobre la influencia del ENOS 
en Norteamérica, ésta presenta una importante 
variación espacial. 
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